Известны ежемесячные продажи австралийского вина в тысячах литров с января 1980 по июль 1995, необходимо построить прогноз на следующие два года.
Попробуем поделить на число дней в месяце: Ряд не стал более регулярным, так что вернёмся к исходным данным.
STL-декомпозиция ряда:
Оптимальное преобразование Бокса-Кокса и результат его применения:
В данном случае преобразование имеет смысл использовать, так как оно хорошо стабилизирует дисперсию. Попробуем округлить параметр и взять \(\lambda=-1\): Результат практически такой же. Далее будем использовать \(\lambda=-1\).
## ETS(A,Ad,A)
##
## Call:
## ets(y = tSeries, lambda = LambdaOpt)
##
## Box-Cox transformation: lambda= -1
##
## Smoothing parameters:
## alpha = 0.0707
## beta = 0.0043
## gamma = 1e-04
## phi = 0.9782
##
## Initial states:
## l = 0.9999
## b = 0
## s=0 0 0 0 0 0
## 0 0 0 0 0 0
##
## sigma: 0
##
## AIC AICc BIC
## -3453.752 -3449.368 -3396.786
Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:
## ME RMSE MAE MPE MAPE MASE
## Training set 51.61557 2189.653 1642.288 -0.3744397 6.478496 0.8362138
## Test set 1232.80660 2127.679 1766.955 3.6014104 6.978830 0.8996913
## ACF1 Theil's U
## Training set -0.08120999 NA
## Test set -0.35511396 0.3246609
Остатки:
Достигаемые уровни значимости критерия Льюнга-Бокса для них:
Остатки не коррелированы
Q-Q plot и гистограмма для остатков:
Распределение имеет выброс слева.
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | отвергается | 1.900071910^{-4} |
| Несмещённость | Уилкоксона | не отвергается | 0.6891151 |
| Стационарность | KPSS | не отвергается | 0.1 |
Применим функцию auto.arima:
## Series: tSeries
## ARIMA(1,1,2)(0,1,2)[12]
## Box Cox transformation: lambda= -1
##
## Coefficients:
## ar1 ma1 ma2 sma1 sma2
## -0.5897 -0.2604 -0.5390 -0.4770 -0.1787
## s.e. 0.6409 0.6149 0.5321 0.0918 0.0866
##
## sigma^2 estimated as 8.281e-08: log likelihood=1780.08
## AIC=-3548.15 AICc=-3547.61 BIC=-3529.63
Предлагается модель ARIMA(1,1,2)(0,1,2)\(_{12}\). Посмотрим на её остатки:
Отрежем первые 13 отсчётов и продолжим анализ:
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | отвергается | 6.367987210^{-5} |
| Несмещённость | Уилкоксона | отвергается | 0.0460378 |
| Стационарность | KPSS | отвергается | 0.0327717 |
Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:
## ME RMSE MAE MPE MAPE MASE
## Training set 543.62191 4285.153 2647.961 2.7109642 11.575578 1.3482784
## Test set -56.43371 1961.410 1517.384 -0.9025861 6.206568 0.7726154
## ACF1 Theil's U
## Training set 0.6188588 NA
## Test set -0.3442693 0.2790998
Исходный ряд нестационарен (p<0.01, критерий KPSS); сделаем сезонное дифференцирование: Ряд всё ещё нестационарен (p<0.01, критерий KPSS). Проведём ещё одно дифференцирование:
Для полученного ряда гипотеза стационарности не отвергается (p>0.1)
Посмотрим на ACF и PACF полученного продифференцированного ряда:
На ACF значимы лаги 1, 12 и 24, на PACF — 1-5. Поищем с помощью auto.arima оптимальную модель полным перебором (stepwise=FALSE) с ограничениями d=1, D=1, max.p=5, max.q=4, max.P=0, max.Q=2, max.order=10:
## Series: tSeries
## ARIMA(1,1,1)(0,1,2)[12]
## Box Cox transformation: lambda= -1
##
## Coefficients:
## ar1 ma1 sma1 sma2
## -0.0201 -0.8693 -0.4843 -0.177
## s.e. 0.0838 0.0415 0.0923 0.087
##
## sigma^2 estimated as 8.229e-08: log likelihood=1780.02
## AIC=-3550.05 AICc=-3549.66 BIC=-3534.61
Была найдена более оптимальная по AICc модель — ARIMA(1,1,1)(0,1,2)\(_{12}\). Посмотрим на её остатки: Отрежем начало ряда остатков и проанализируем их:
Достигаемые уровни значимости критерия Льюнга-Бокса для остатков:
Q-Q plot и гистограмма:
| Гипотеза | Критерий | Результат проверки | Достигаемый уровень значимости |
|---|---|---|---|
| Нормальность | Шапиро-Уилка | отвергается | 4.503317210^{-5} |
| Несмещённость | Уилкоксона | не отвергается | 0.055308 |
| Стационарность | KPSS | отвергается | 0.0382442 |
Настроив выбранную модель на обучающей выборке, посчитаем её качество на тестовой:
## ME RMSE MAE MPE MAPE MASE
## Training set 535.674930 4312.632 2699.617 2.6262660 11.786210 1.3745806
## Test set 4.690032 1936.939 1490.400 -0.7019034 6.133156 0.7588758
## ACF1 Theil's U
## Training set 0.6172789 NA
## Test set -0.3603302 0.2778012
Сравним остатки двух версий аримы, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:
##
## Diebold-Mariano Test
##
## data: resres.auto
## DM = 0.21196, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.8324
## alternative hypothesis: two.sided
Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов.
В целом подобранная вручную модель проще, её остатки лучше, а ошибка на тесте меньше, так что остановимся на модели, подобранной вручную.
Сравним остатки лучших моделей ARIMA и ETS, одинаково обрезав их начало так, чтобы у обоих методов они были правильно определены:
##
## Diebold-Mariano Test
##
## data: res.autores.ets
## DM = 1.4329, Forecast horizon = 1, Loss function power = 2,
## p-value = 0.1538
## alternative hypothesis: two.sided
Критерий Диболда-Мариано не обнаруживает значимого различия между качеством прогнозов, но у ARIMA ошибка на тесте и AIC меньше, так что выберем её.
Поскольку остатки ненормальны, предсказательные интервалы строятся бутстрепом
## Warning in InvBoxCox(pred$pred, lambda, biasadj, var(residuals(object), :
## biasadj information not found, defaulting to FALSE.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Aug 1994 29266.72 2488.820 NA 1676.706 NA
## Sep 1994 24226.65 2432.219 NA 1647.597 NA
## Oct 1994 27958.44 2447.141 NA 1650.091 NA
## Nov 1994 32727.79 2460.692 NA 1651.952 NA
## Dec 1994 37983.85 2468.522 NA 1651.224 NA
## Jan 1995 15250.22 2235.830 NA 1540.085 NA
## Feb 1995 21981.02 2325.168 NA 1578.128 NA
## Mar 1995 24016.60 2330.496 NA 1576.790 NA
## Apr 1995 25868.30 2331.197 NA 1573.366 NA
## May 1995 24155.14 2301.301 NA 1556.054 NA
## Jun 1995 25890.01 2301.108 NA 1552.372 NA
## Jul 1995 30204.94 2315.639 NA 1555.389 NA
## Aug 1995 28523.51 2004.960 NA 1343.665 NA
## Sep 1995 24887.50 1965.200 NA 1321.083 NA
## Oct 1995 27791.14 1960.512 NA 1313.995 NA
## Nov 1995 32830.91 1961.203 NA 1309.437 NA
## Dec 1995 38462.01 1958.033 NA 1303.254 NA
## Jan 1996 15994.05 1810.062 NA 1231.788 NA
## Feb 1996 21658.53 1847.446 NA 1244.730 NA
## Mar 1996 24290.08 1847.047 NA 1240.366 NA
## Apr 1996 25943.91 1838.750 NA 1232.530 NA
## May 1996 24453.66 1814.179 NA 1217.492 NA
## Jun 1996 25509.52 1803.462 NA 1208.801 NA
## Jul 1996 30508.44 1808.362 NA 1207.192 NA